#include <TRandom3.h>
#include <numeric>
const float Me_MUL_C2 = 0.511*1e6;//eV

int ceilToInt(const float x) 
{
  return static_cast<int>(ceil(x));
}

float calculateAverage(const std::vector<float>& data) {
    if (data.empty()) {
        return 0.0f; 
    }

    float sum = std::accumulate(data.begin(), data.end(), 0.0f);
    return sum / data.size();
}

TTree* getEventTree(string fileName, string treeName)
{   
    TFile* inputFile = new TFile(fileName.c_str());
    TTree* eventTree = (TTree*)inputFile->Get(treeName.c_str());
    cout << "--> Root file loaded : " << fileName << endl;
    int num_of_entries = eventTree->GetEntries();
    cout << "the number of total entries is : " << num_of_entries << endl;

    return eventTree;
}

float dot_product(float a[], float b[]) 
{
    float result = 0;
    for (int i = 0; i < 3; i++) {
        result += a[i] * b[i];
    }
    return result;
}

float module_length(float a[3]) 
{
    float result = 0;
    for (int i=0; i < 3; i++)
    {
        result += pow(a[i], 2);
    }
    return sqrt(result);
}

void normalize_vector(float a[3])
{
    float mag_a = module_length(a);
    for (int i=0; i < 3; i++)
    {
        a[i] = a[i]/mag_a;
    }
}

void mul_vector(float a[3], float length)
{
    for (int i=0; i < 3; i++)
    {
        a[i] = a[i] * length;
    }
}

float angle_between_vectors(float a[3], float b[3]) 
{
    float dot = dot_product(a, b);
    float mag_a = module_length(a);
    float mag_b = module_length(b);
    float cos_theta = dot / (mag_a * mag_b);
    float theta = acos(cos_theta) * 180 / M_PI;
    // cout << " theta_r: " << theta << "  ";
    return theta;
}

TRandom3 fRandom;
float smearing(float energy_actual, float energyRef, float resolution_at_Eref)
{   
    if (energy_actual > 0)
    {
        float eRes = resolution_at_Eref * sqrt(energyRef / energy_actual) / 2.35 / 100.0;
        float gain = fRandom.Gaus(1.0, eRes);
        float energy_smearing = energy_actual *gain;
        return energy_smearing;
    }
    return 0;
}


float camp_calculate(float energy_elec, float energy_scatter)
{   
    energy_elec = energy_elec*1e3;
    energy_scatter = energy_scatter*1e3;
    float cos_em = 1 + Me_MUL_C2/(energy_elec + energy_scatter) - Me_MUL_C2/energy_scatter;
    // cout << " " << Me_MUL_C2/(energy_elec + energy_scatter) << " " << Me_MUL_C2/energy_scatter << " ";
    cos_em = (cos_em > 1.0) ? 1.0 : ((cos_em < -1.0) ? -1.0 : cos_em);
    // cout << "cos_em: " << cos_em << " " << endl;
    float theta = acos(cos_em) * 180 / M_PI;
    // cout << "theta_em: " << theta << "  ";
    return theta;
}


float calculate_angle_error(float e_vector[3], float electron_energy, float gamma_vector[3], float gamma_energy, float vector_actual[3])
{
    // cout << "electron_energy: " << electron_energy ;
    electron_energy = sqrt(pow(electron_energy, 2) + electron_energy*1.022);
    // cout << " electron_p: " << electron_energy << endl;
    normalize_vector(e_vector);
    cout << "e_vector: (" << e_vector[0] << ", " << e_vector[1] << ", "<< e_vector[2] << ")" << endl;
    mul_vector(e_vector, electron_energy);
    normalize_vector(gamma_vector);
    mul_vector(gamma_vector, gamma_energy);
    float vector_incidence_calculate[3];
    for (int i=0; i < 3; i++)
    {
        vector_incidence_calculate[i] = round((e_vector[i] + gamma_vector[i])*100)/100;
    }
    normalize_vector(vector_incidence_calculate);
    // float angle_error = round(angle_between_vectors(vector_incidence_calculate, vector_actual)*10)/10;
    // if (vector_incidence_calculate[0] < 0) angle_error = angle_error*(-1);
    cout << "vector_incidence_calculate:("  << vector_incidence_calculate[0] << ", "
    << vector_incidence_calculate[1] << ", " << vector_incidence_calculate[2] << ")" << endl;
    // << "angle_error: " << angle_error << endl;
   
    return (vector_incidence_calculate[0] - vector_actual[0]);
    // return angle_error;
}

float calculate_angle_error_y(float e_vector[3], float electron_energy, float gamma_vector[3], float gamma_energy, float vector_actual[3])
{
    // cout << "electron_energy: " << electron_energy ;
    electron_energy = sqrt(pow(electron_energy, 2) + electron_energy*1.022);
    // cout << " electron_p: " << electron_energy << endl;
    normalize_vector(e_vector);
    cout << "e_vector: (" << e_vector[0] << ", " << e_vector[1] << ", "<< e_vector[2] << ")" << endl;
    mul_vector(e_vector, electron_energy);
    normalize_vector(gamma_vector);
    mul_vector(gamma_vector, gamma_energy);
    float vector_incidence_calculate[3];
    for (int i=0; i < 3; i++)
    {
        vector_incidence_calculate[i] = round((e_vector[i] + gamma_vector[i])*100)/100;
    }
    normalize_vector(vector_incidence_calculate);
    // float angle_error = round(angle_between_vectors(vector_incidence_calculate, vector_actual)*10)/10;
    // if (vector_incidence_calculate[0] < 0) angle_error = angle_error*(-1);
    // cout << "vector_incidence_calculate:("  << vector_incidence_calculate[0] << ", "
    // << vector_incidence_calculate[1] << ", " << vector_incidence_calculate[2] << ")" << endl;
    // << "angle_error: " << angle_error << endl;
   
    return (vector_incidence_calculate[1] - vector_actual[1]);
    // return angle_error;
}


int getMainTrackID(TRestG4Event* g4Event)
{   
    int fNumOfTrack = g4Event->GetNumberOfTracks();
    // cout << "fNumOfTrack: " << fNumOfTrack << endl;
    if (fNumOfTrack==0) return 0;
    int fTrackID = 0;
    int fMainfTrackID = 0;
    TRestG4Track* fTrack = new TRestG4Track();
    int max_hits = 0;
    int fNumOfHits = 0;
    for (int i=0; i < fNumOfTrack; i++)
    {
        fTrackID = i;
        fTrack = g4Event->GetTrack(i);
        fNumOfHits = fTrack->GetNumberOfHits();
        // cout << "track_" << fTrackID << ": " << fNumOfHits << endl;
        if (fNumOfHits >= max_hits)
        {
            fMainfTrackID = fTrackID;
            max_hits = fNumOfHits;
        }
    }
    // cout << "fMainfTrackID: " << fMainfTrackID << endl;
    return fMainfTrackID;
}

void travelAllHits(TRestG4Event* g4Event)
{
    int fNumOfTrack = g4Event->GetNumberOfTracks();
    TRestG4Track* fTrack = new TRestG4Track();
    TRestG4Hits* fHits = new TRestG4Hits();
    int fTrackID = 0;
    int fNumOfHits = 0;
    for (int i=0; i < fNumOfTrack; i++)
    {
        fTrackID = i;
        fTrack = g4Event->GetTrack(i);
        fNumOfHits = fTrack->GetNumberOfHits();
        fHits = fTrack->GetHits();
        for (int j=0; j < fNumOfHits; j++)
        {
            cout << "Track_" << i << ", Hits_" << j << " :(" << fHits->GetX(j) << ", " << fHits->GetY(j) << ", " << fHits->GetZ(j) << ") " 
            << fHits->GetEnergy(j) << endl;
        }
        cout << endl;
    }
}

float getTrackEnergy(TRestG4Track* fTrack)
{
    int fNumOfHits = fTrack->GetNumberOfHits();
    TRestG4Hits* fHits = new TRestG4Hits();
    fHits = fTrack->GetHits();
    float fTrackTotalEnergy = 0.;
    for (int j=0; j < fNumOfHits; j++)
    {
        fTrackTotalEnergy += fHits->GetEnergy(j);
    }
    return fTrackTotalEnergy;
}

float calculate_distance_two_hits(float hit0[3], float hit1[3])
{
    float pow_distance = 0.;
    for(int i=0; i<3; i++)
    {
        pow_distance += pow((hit0[i]-hit1[i]), 2);
    }
    return sqrt(pow_distance);
}
